import cv2 as cv
import matplotlib.pyplot as plt
import numpy as np
import copy

img = cv.pyrDown(cv.imread(r'img/Lena.jpg'))
gray = cv.cvtColor(img, cv.COLOR_BGR2GRAY)

print('第一题')
img_blur = cv.blur(img, (5, 5))
img_g_blur = cv.GaussianBlur(img, (5, 5), 5)


def imshow1():
    cv.imshow('Lena', img)
    cv.imshow('Blur', img_blur)
    cv.imshow('Gaussian Blur', img_g_blur)
    if cv.waitKey(0) == ord('c'):
        cv.destroyAllWindows()


# imshow1()


print('第二题')
x = cv.Sobel(gray, cv.CV_16S, 1, 0)
y = cv.Sobel(gray, cv.CV_16S, 0, 1)
absX = cv.convertScaleAbs(x)  # 转回uint8
absY = cv.convertScaleAbs(y)
sobel = cv.addWeighted(absX, 0.5, absY, 0.5, 0)
canny = cv.Canny(gray, 100, 150)


def imshow2():
    cv.imshow("Sobel", sobel)
    cv.imshow("Canny", canny)
    if cv.waitKey(0) == ord('c'):
        cv.destroyAllWindows()


# imshow2()


print('第三题')
rice = cv.pyrDown(cv.imread(r'img/rice.jpg'))
rice_gray = cv.cvtColor(rice, cv.COLOR_BGR2GRAY)
hist_rice = cv.calcHist([rice_gray], [0], None, [256], [0, 256])
rice_gray_t = cv.threshold(rice_gray, 0, 255, cv.THRESH_OTSU)


def imshow3():
    plt.subplot(121), plt.imshow(rice_gray, cmap='gray')
    plt.subplot(122), plt.plot(hist_rice)
    plt.show()
    cv.imshow("threshold_otsu", rice_gray_t[1])
    if cv.waitKey(0) == ord('c'):
        cv.destroyAllWindows()


# imshow3()

print('第四题')

''' 
https://blog.csdn.net/u011276025/article/details/89790190
opencv 一种不均匀光照的补偿方法
'''


def unevenLightCompensate(img, blockSize):
    gray = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
    average = np.mean(gray)

    rows_new = int(np.ceil(gray.shape[0] / blockSize))
    cols_new = int(np.ceil(gray.shape[1] / blockSize))

    blockImage = np.zeros((rows_new, cols_new), dtype=np.float32)
    for r in range(rows_new):
        for c in range(cols_new):
            rowmin = r * blockSize
            rowmax = (r + 1) * blockSize
            if (rowmax > gray.shape[0]):
                rowmax = gray.shape[0]
            colmin = c * blockSize
            colmax = (c + 1) * blockSize
            if (colmax > gray.shape[1]):
                colmax = gray.shape[1]

            imageROI = gray[rowmin:rowmax, colmin:colmax]
            temaver = np.mean(imageROI)
            blockImage[r, c] = temaver

    blockImage = blockImage - average
    blockImage2 = cv.resize(blockImage, (gray.shape[1], gray.shape[0]), interpolation=cv.INTER_CUBIC)
    gray2 = gray.astype(np.float32)
    dst = gray2 - blockImage2
    dst = dst.astype(np.uint8)
    # dst = cv.GaussianBlur(dst, (3, 3), 0)  # 去掉原方法中的高斯模糊避免影响面积计算
    dst = cv.cvtColor(dst, cv.COLOR_GRAY2BGR)

    return dst


'''
为避免亮度不均匀对阈值分割产生影响，先对图像进行亮度均匀化，然后用腐蚀分离黏连的米粒再计算数量。
由于腐蚀后米粒面积和长度会减小，因此非黏连的米粒用回腐蚀前的图像进行面积和长度计算。
'''
rice_compensated = unevenLightCompensate(rice, 50)
rice_com_gray = cv.cvtColor(rice_compensated, cv.COLOR_BGR2GRAY)
rice_threshold = cv.threshold(rice_com_gray, 0, 255, cv.THRESH_OTSU)[1]
kernel = cv.getStructuringElement(cv.MORPH_RECT, (3, 3))  # 腐蚀卷积核
rice_erode = cv.erode(rice_threshold, kernel)
element = cv.getStructuringElement(cv.MORPH_CROSS, (5, 5))
rice_a = cv.morphologyEx(rice_erode, cv.MORPH_OPEN, element)
rice_b = cv.morphologyEx(rice_threshold, cv.MORPH_OPEN, element)
seg = copy.deepcopy(rice_a)
seg2 = copy.deepcopy(rice_b)
binary, contours, hierarchy = cv.findContours(seg, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE)
# contours, hierarchy = cv.findContours(seg, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE)
binary2, contours2, hierarchy2 = cv.findContours(seg2, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE)
# contours2, hierarchy2 = cv.findContours(seg2, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE)
count = 0
size_list = []
area_list = []
length_list = []
info_list = []

'''由于腐蚀后米粒面积和长度会减小，因此非黏连的米粒用回腐蚀前的图像进行面积和长度计算。'''
for i in range(len(contours2), 0, -1):
    c = contours2[i - 1]
    area = cv.contourArea(c)
    if area < 50:
        continue
    # count = count + 1
    x, y, w, h = cv.boundingRect(c)
    length = np.sqrt(w ** 2 + h ** 2)
    size_list.append((x, y, area, length))
    # area_list.append(area)
    # length_list.append(length)
    # info_list.append((area, length))

for i in range(len(contours), 0, -1):
    c = contours[i - 1]
    area = cv.contourArea(c)
    if area < 50:
        continue
    count = count + 1

    x, y, w, h = cv.boundingRect(c)
    length = np.sqrt(w ** 2 + h ** 2)
    for size in size_list:
        if abs(x - size[0]) < 5 and abs(y - size[1]) < 5:
            if area * 2 > size[2]:
                area = size[2]
                length = size[3]
    area_list.append(area)
    length_list.append(length)
    info_list.append((area, length))
    print("blob", count, ", area:", area, "length:", length)
    cv.rectangle(rice, (x, y), (x + w, y + h), (0, 0, 0xff), 1)
    x = x if x < 12 else x - 8
    y = y if y > 12 else y + 8
    cv.putText(rice, str(count), (x, y), cv.FONT_HERSHEY_PLAIN, 0.8, (0, 0xff, 0))

area_mean = np.mean(area_list)
length_mean = np.mean(length_list)
area_std = np.std(area_list)
length_std = np.std(length_list)
print("米粒数量：", count)
print("面积均值：", np.round(area_mean, 2))
print("长度均值：", np.round(length_mean, 2))
print("面积方差：", np.round(np.var(area_list), 2))
print("长度方差：", np.round(np.var(length_list), 2))
valid_area_min = round(area_mean - 3 * area_std, 2)
valid_area_max = round(area_mean + 3 * area_std, 2)
valid_length_min = round(length_mean - 3 * length_std, 2)
valid_length_max = round(length_mean + 3 * length_std, 2)
print("面积有效范围：", valid_area_min, "-", valid_area_max)
print("长度有效范围：", valid_length_min, "-", valid_length_max)
valid_count = 0

for info in info_list:
    if valid_area_max > info[0] > valid_area_min and valid_length_max > info[1] > valid_length_min:
        valid_count = valid_count + 1

print("落在 3 sigma 内米粒数量：", valid_count)


def imshow4():
    cv.imshow("rice_threshold", rice_a)
    cv.imshow("rice_threshold_erode", rice_b)
    cv.imshow("rice", rice)
    if cv.waitKey(0) == ord('c'):
        cv.destroyAllWindows()


# imshow4()


print('第五题')


def my_detect(name, method, image, image_gray):
    kp = method.detect(image_gray, None)
    out = None
    out = cv.drawKeypoints(image, kp, out)
    cv.putText(out, name, (20, 20), cv.FONT_HERSHEY_PLAIN, 2, (0xff, 0xff, 0))
    return out


chess = cv.imread(r'img/left03.jpg')
road = cv.pyrDown(cv.imread(r'img/road.jpg'))
chess_gray = cv.cvtColor(chess, cv.COLOR_BGR2GRAY)
road_gray = cv.cvtColor(road, cv.COLOR_BGR2GRAY)
# SIFT
sift = cv.xfeatures2d.SIFT_create()
chess_sift = my_detect("SIFT", sift, chess, chess_gray)
road_sift = my_detect("SIFT", sift, road, road_gray)
# FAST
fast = cv.FastFeatureDetector_create(threshold=40, nonmaxSuppression=True, type=cv.FAST_FEATURE_DETECTOR_TYPE_9_16)
chess_fast = my_detect("FAST", fast, chess, chess_gray)
road_fast = my_detect("FAST", fast, road, road_gray)
# ORB
orb = cv.ORB_create()
chess_orb = my_detect("ORB", orb, chess, chess_gray)
road_orb = my_detect("ORB", orb, road, road_gray)

chess_result = np.concatenate([chess_sift, chess_fast, chess_orb], axis=1)
road_result = np.concatenate([road_sift, road_fast, road_orb], axis=1)
chess_result = cv.resize(chess_result, (int(chess_result.shape[1] * 0.8), int(chess_result.shape[0] * 0.8)))
road_result = cv.resize(road_result, (int(road_result.shape[1] * 0.8), int(road_result.shape[0] * 0.8)))
# result = np.concatenate([chess_result, road_result], axis=1)


# Harris
chess_gray_float = np.float32(chess_gray)
dst = cv.cornerHarris(chess_gray_float, 2, 3, 0.04)
dst = cv.dilate(dst, None)
chess[dst > 0.3 * dst.max()] = [0, 0, 255]


def imshow5():
    cv.imshow("chess_Harris", chess)
    cv.imshow("chess_result", chess_result)
    cv.imshow("road_result", road_result)


imshow5()

if cv.waitKey(0) == ord('q'):
    cv.destroyAllWindows()
